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RESUMEN 


Presentamos simulaciones numericas de la interaccion entre un “pulso” bipo¬ 
lar colimado expulsado por una estrella y el viento continuo debido a una estrella 
compafiera. Se muestra la variacion de las propiedades de los mapas de intensidad 
en Ha al variar los valores de algunos de los parametros de entrada. Mostramos 
que la asimetria (en el tamafio y la intensidad) entre los dos lobulos de la nebulosa 
proto-planetaria OH 231.8-h4.2 se puede explicar en el escenario propuesto si el 
viento de la estrella compafiera muestra una fuerte dependencia con la latitud. 

ABSTRACT 

We present numerical simulations of the interaction between a collimated, 
bipolar “pulse” ejected from a star and a continuous wind ejected from a stellar 
companion. We explore the characteristics of the predicted Ha intensity maps 
by varying selected input parameters. We find that the asymmetry (in size and 
strength) between the two lobes of the proto-planetary nebula OH 231.8-1-4.2 is 
reproduced in this scenario if the wind ejected by the companion star has a strong 
latitude dependence. 

Key Words: ISM: KINEMATICS AND DYNAMICS - PLANETARY 
NEBULA: INDIVIDUAL (OH231.8+4.2) 


1. INTRODUCTION 

Most proto-planetary Nebulae (PPNe) show 
bipolar or mutipolar lobes, highly collimated out¬ 
flows or jets and other complex structures. A partic¬ 
ularly well studied PPN candidate is the nebula as¬ 
sociated with the OH 231.8+4.2 source. This is a re¬ 
markable example of a bipolar PPNe with collimated 
outflows. In Ha and forbidden emission line maps, 
this object presents a clear axial symmetry with the 
existence of shocked material in a wide, bipolar bub¬ 
ble with an obscuring ridge in between the two lobes 
(Bujarrabal et al. 2002). Long-slit spectra of this 
object show the presence of Herbig-Haro like knots 
at the end of the bipolar lobes, which probably trace 
the interaction between the fast, collimated molecu¬ 
lar jet and the ambient material (Sanchez Contreras 
et al. 2000a). Contrary to what happens in most 
PPNe, the two optical lobes have different sizes; the 
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Southern lobe is weaker but about two times more 
extended than the Northern lobe (Reipurth 1987; 
Sanchez Contreras et al. 2000a). 

In the optical and NIR continuum, the nebula is 
narrow and elongated along its symmetry axis (Al¬ 
colea et al. 2001; Meakin et al. 2003). The same nar¬ 
row component is observed in CO and other molecu¬ 
lar emission lines (Alcolea et al. 2001; Sanchez Con¬ 
treras, Bujarrabal & Alcolea 1997; Sanchez Contr¬ 
eras et al. 2000b). The kinematical properties of 
this molecular jet are remarkable; its velocity in¬ 
creases linearly with distance from the star following 
a Hubble-like velocity law, which is consistent with 
a sudden acceleration event (Alcolea et al. 2001; Al¬ 
colea 2004). The Southern molecular jet (as seen 
in CO emission) is clearly more extended than the 
Northern jet and shows the highest flow velocity. 
The deprojected velocities show values of up to 430 
km s“^ in the Southern jet, and up to 210 km s“^ in 
the Northern jet (Alcolea et al. 2001). The length 
along the axis of the OH 231.8+4.2 nebula is ^50", 
which at a distance of 1.5 kpc (Kastner et al. 1992) 
implies a length >10^^ cm. 

The central source of OH 231.8+4.2 is not de- 
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tected at optical wavelengths because of the dusty 
envelope, but the reflected light has a M9 III spec¬ 
tral type, with a blue excess that suggests the pres¬ 
ence of a warmer companion (or a peculiar photo- 
spheric structure and/or an accretion disk; Cohen 
1981). This star is a typical Mira variable at the end 
of the Asymptotic Giant Branch (AGB) identified 
with QX Pup (Kastner et ah 1998). Recently, OH 
231.8+4.2 has been mapped at high spatial resolu¬ 
tion in SiO and H 2 O maser emission (Sanchez Con¬ 
treras et al. 2002; Gomez & Rodriguez 2001). The 
position of the SiO and H 2 O masers (which are trac¬ 
ing the position of QX Pup), are clearly offset from 
the axis of the bipolar outflow (as first pointed out 
by Gomez & Rodriguez 2001). 

The QX Pup star shows the variability and maser 
emission properties typical of a Mira star (or of an 
OH/IR star) at the AGB. The contradiction between 
the apparent evolutionary stage of the nebula (i.e., 
PPN) and that of the QX Pup (AGB) star has been 
discussed by many authors (see, e.g., Alcolea et al. 

2001) . Based on the large displacement (> 1000 AU) 
between the position of the SiO maser and the bipo¬ 
lar axis of the nebula, Alcolea (2004) has proposed 
that QX Pup may not be directly related to the bipo¬ 
lar nebula. As first pointed out by Alcolea (2004), if 
the QX Pup star were located in the Northern lobe 
the observed differences between the two lobes possi¬ 
bly could be explained as a result of the interaction 
between the northern lobe and the wind from QX 
Pup. 

In this paper, we present 3D numerical simula¬ 
tions of the interaction between a collimated, bipo¬ 
lar “pulse” ejected from a star and a continous wind 
ejected from a stellar companion. These simula¬ 
tions are a first attempt to reproduce the asym¬ 
metry between the two lobes of the bipolar PPN 
OH 238.1+4.2 in terms of a bipolar outflow/stellar 
wind interaction model. Erom the numerical simula¬ 
tions, we obtain predictions of Ra maps that can be 
compared with the optical images of OH 238.1+4.2 
(Sanchez Gontreras et al. 2000a; Buiarrabal et al. 

2002 ) . 

The paper is organized as follows. In § 2 we 
describe the parameters and the numerical method 
which have been used. In § 3 we present the results 
of the numerical simulations. In § 4 we compare our 
simulations with the observations of the PPN OH 
231.8+4.2. 

2. THE NUMERIGAL MODELS 
2.1. General eonsiderations 

We compute a series of five models of the inter¬ 
action between a collimated, bipolar “pulse” ejected 


from a star and a continuous wind ejected from a 
stellar companion. The models are computed in a 
3D cartesian coordinate system, with its origin cen¬ 
tered in the middle of the computational domain. 
The domain has a physical extent of 1.5 x 10^^ cm 
along the z-axis, and of 7.5 x 10^^ cm along the x- 
and ^-axes. Outflow conditions are applied in all of 
the domain boundaries. 

The computations are carried out with the 
“yguazu-a” code, which integrates the gasdynamic 
equations in a binary, adaptive computational grid. 
This code is described in detail by Raga, Villagran- 
Muniz & Navarro-Gonzalez (2000). In the version of 
the code that we have employed, a single rate equa¬ 
tion for neutral hydrogen is integrated (together with 
the 3D gasdynamic equations), and a simplified cool¬ 
ing function is calculated as a function of the density, 
ionization fraction and temperature (as described by 
Raga et al. 2000). This cooling function is then 
added as a sink term in the energy equation. 

Eor our simulations, we have used a 5-level binary 
grid with a maximum resolution of 5.86 x 10^^ cm 
along the three axes. The domain was initialized 
with a wind solution (§2.2) and a “bipolar pulse” 
(§2.3), and the simulations were carried out until 
the shell pushed out by the “pulse” reaches one of 
the boundaries of the computational domain. The 
initial flow configuration is shown in Eigure^ 

2.2. The wind 

We have initialized the computational domain 
with a wind with a latitude-dependent density of the 
parametrized form : 


P = 


47rr^ 


f{0): 


( 1 ) 


where r is the spherical radius measured from the po¬ 
sition of the wind source. The wind has an outward 
directed velocity (independent of r) given by : 


V = 


Vw 

/W’ 


( 2 ) 


with Vw = 40 km s“^ for all of the computed models, 
and with different values of Mw for the successive 
models (see Table 1). The anysotropy function f{0) 
has the parametrized form : 


/W=e-(e-l)cos^^, (3) 


where 0 is the angle measured from the polar axis of 
the wind (which we take to be parallel to the z-axis 
of the computational domain), ^ is a constant which 
is equal to the equator-to-pole density ratio, and p is 
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Fig. 1. Schematic diagram showing the configuration 
used for the numerical simulations (see the text). A 
collimated wind “pulse” is initialized within a sphere of 
radius Tc (centered at the origin of the coordinate sys¬ 
tem) and the continuous wind from the stellar compan¬ 
ion is imposed (at all times) within a sphere of radius Tw , 
which is centered at the point {xw,0, Zw)- In the diagram, 
the projections of these two spheres on the xz-plane are 
shown. The axis of the “pulse” (as well as the axis of the 
wind in the simulations with a latitude-dependent wind) 
is parallel to the z-axis. 

a constant which determines the degree of flattening 
towards the equator of the density stratification. As 
can be seen from Equation (|3|), a wind with ^ = 1 is 
isotropic. The choices of ^ and p for the computed 
models are given in Table 1. 

The source of the wind is located at a position 
away from the origin of the coordinate 
system (which is at the center of the computational 
domain, see EigureHand §2.1). The values of and 
Zw chosen for the computed models are given in Ta¬ 
ble 1. The whole of the computational domain (with 
the exception of the region occupied by the “bipolar 
pulse”, see §2.3) is initialized with the wind solution, 
with the density and velocity given by Equations m 
ini. Also, the wind solution is imposed at all times 
within a sphere of radius = 2.5 x 10 ^^ cm, cen¬ 
tered on the {xw^{)^Zw) source position (see Eigure 

ID. 

We impose a uniform, = 10 ^ K within the 
sphere of radius , and a temperature 

T = (4) 


TABLE 1 

COLLIMATED PULSE/WIND INTERACTION 
MODELS 


Model 

[ 10 “®Moyr“i] 

rf^ b 

[IQi" 

cm] 

r 


[yr] 

Ml 

1 

- 1.0 

1.5 

1 

— 

1100 

M2 

5 

- 1.0 

1.5 

1 

— 

1500 

M3 

3 

- 1.0 

1.5 

20 

0.5 

1500 

M4 

10 

-3.0 

1.5 

20 

0.5 

2300 

M5 

10 

-3.0 

1.5 

50 

0.5 

2700 


^Mass loss of the wind from the stellar companion 
^Position on the xz-plane of the stellar wind source 
^Asymmetry parameters of the wind 

*^Time at which the expanding shell starts to leave the 
computational grid (the frames of figures 2-6 correspond 
to these times) 

outside of this sphere, where and are con¬ 
stants defined above, 7 = 5/3 is the specific heat 
ratio and r is the spherical radius measured from 
the position of the wind source. This temperature 
stratification (corresponding to a constant velocity, 
adiabatic wind) is of course only used in the r > 
region as an initial condition. Einally, the ejected 
wind is assumed to have fully neutral H, and to have 
a small, seed ionization fraction coming from singly 
ionized C. 

2.3. The bipolar pulse 

As an initial condition, we impose a “bipolar 
pulse” of material ejected from a source located at 
the origin (x, z) = ( 0 , 0 , 0 ) of the coordinate sys¬ 
tem (see Eigure nj. This bipolar “pulse” is imposed 
as an initial condition within a sphere of radius Vc 
(centered on the origin, see Eigure ^), and has a 
“hubble-like velocity law”, radially directed velocity 
of the form 



where R is the spherical radius (measured from the 
position of the bipolar outflow source), and a density 



p = Pc g{v), 

(6) 

where 

Vc sin Oc 

(7) 


RsinO ’ 

and 

g{T]) = min[l,7?]. 

( 8 ) 


Equations iIHlHIi give a cylindrically stratified density 
distribution, with a core of cylindrical radius Vc sin Oc 
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of constant density pc, and densities that fall mono- 
tonically for larger cylindrical radii. This kind of 
density stratification for the jet-like flow is inspired 
in the asymptotic MHD wind solution of Shu et ah 
(1995). 

In all of our models we have considered Vc = 
5 X 10^^ cm, Uc = Pc /= 10^ cm“^ (where 
rriH is the mass of the hydrogen atom) for the cen¬ 
tral density of the “collimated pulse” (equation |n|) , 
Oc = 5° (equation [3) and Vc = 200 km s“^. Eor 
the “pulse” material (within the sphere of radius Vc, 
see above and figure Q) we have assumed a uniform 
temperature of 100 K, and that the gas is neutral 
(except for a seed ionization fraction coming from 
singly ionized C). 

The bipolar “pulse” is then “released” at t = 0 
(i. e., the material within the sphere of radius Vc is 
allowed to evolve freely as a function of time), and 
the gasdynamic equations are integrated in time to 
follow the interaction of the “pulse” with the stellar 
wind (described in §2.2). The results of the numer¬ 
ical integrations are described in the following sec¬ 
tion. 


3. MODEL RESULTS 

We have then computed 5 models (models MI¬ 
MS), with identical “bipolar pulses” (see §2.3) inter¬ 
acting with winds of different parameters (coming 
from the binary companion of the source of the bipo¬ 
lar pulse, see Eigure^. We have tried both isotropic 
winds (models Ml and M2) and latitude dependent 
winds (models M3-M5, see table 1), as described in 
§ 2 . 2 . 

In Eigures EE we show some of the results ob¬ 
tained from models M1-M5 for the times tf (see Ta¬ 
ble 1) at which the shell which is pushed out by 
the bipolar “pulse” reaches one of the boundaries 
of the computational domain. These Eigures show 
the density stratifications on the xz-plane (which in¬ 
cludes the positions of both the wind and the “bipo¬ 
lar pulse” sources, see EigureOJ and on the ^ 2 ;-plane. 
These Eigures also show the Ha maps obtained by 
integrating the Ha emission coefficient (including the 
radiative recombination cascade and the n = 1 ^ 3 
collisional excitations) along the y- and x-axes. 

The models with spherically symmetric winds 
(models Ml and M2, Table 1) produce highly asym¬ 
metric xz-plane density stratifications and xz-Ha 
emission maps (left hand side frames of EiguresEland 
El). The upper lobe of the shell (which is pushed out 
by the bipolar pulse) shows a major “indentation” 
formed by the interaction with the dense region at 
the base of the stellar wind. The ^z-plane cuts and 



X [cm] 


y [cm] 


Fig. 2. Results obtained from model Ml for a t/ = 
1100 yr (see the text and table 1). The top frames show 
the density stratifications on the xz- (left) and yz-plane 
(right) with a logarithmic greyscale (given in g cm“^ by 
the bar on the right). The bottom frames show the Ho 
maps obtained by integrating the Ho emission coefficient 
along the y- (left) and x-axis (right) with a logarithmic 
greyscale (given in erg cm“^ s“^ sterad”^ by the bar on 
the right). 

Ho maps (right hand side frames of Eigures El and 
E|) show an upper lobe which is narrower and shorter 
than the bottom lobe, with a larger asymmetry for 
the model with larger mass loss rate (model M2, see 
table 1). 

In model M3, we have a wind mass loss rate in¬ 
termediate between the ones of models Ml and M2 
(see Table 1), but we have considered a latitude- 
dependent wind (with ^ = 20 and p = 0.5, see equa- 
tiouEl)- In model M3 (see Eigure^, the bottom lobe 
of the outflow is somewhat narrower, but otherwise 
similar to the corresponding lobes of models Ml and 
M2 (see Eigures El and E|)- The upper lobe of model 
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X [cm] y [cm] 

Fig. 3. Results obtained from model M2 for a t/ = 
1500 yr (see the text and table 1). A description of the 
four frames is given in the caption of Figure 1. 

M3, however, is strongly reduced in size, leading to a 
large asymmetry between the top and bottom lobes 
of the expanding shell. 

In models M4 and M5, we have moved the stel¬ 
lar wind source away from the axis of the “bipo¬ 
lar pulse” ejection (see Table 1 and Figures 0 and 
|H|). These models have a high, = lO“^M 0 yr“^ 
mass loss rate and two different degrees of asymme¬ 
try (^ = 20 and 50 for models M4 and M5, respec¬ 
tively). These models produce top-to-bottom lobe 
asymmetries which are qualitatively similar to the 
ones obtained from model M3 (see above and Figure 

El). 

As can be seen from the left hand side plots of 
Figures [5| and El models M4 and M5 have the in¬ 
teresting property that the side-to-side asymmetries 
of the top lobe (when seen in the xz-plane density 
cut or Ha map) are much less important than the 



X [cm] y [cm] 

Fig. 4. Results obtained from model M3 for a t/ = 
1500 yr (see the text and table 1). A description of the 
four frames is given in the caption of Figure 1. 

ones found in models M1-M3. This result indicates 
that a wind with a strong latitude dependence from 
a source placed well away from the axis of the “bipo¬ 
lar pulse” produces a strong asymmetry between the 
two lobes of the expanding shell, but does not intro¬ 
duce strong deviations from the initial axisymmetry 
of the bipolar pulse. 

Models M4 and M5 show a double peak struc¬ 
ture at the bottom lobe (see Figures 0 and EJ, which 
are reminiscent of the results obtained from numer¬ 
ical simulations of stellar jets for a high jet to envi¬ 
ronment density ratio (e.g., Raga 1988; Downes & 
Cabrit 2003). 

4. DISCUSSION 

The results of our 3D numerical simulations of 
the interaction between a bipolar “pulse” and a con- 
tinous wind ejected by a companion star indicate 
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Fig. 5. Results obtained from model M4 for a t/ = 
2300 yr (see the text and table 1). A description of the 
four frames is given in the caption of Figure 1. 


Fig. 6. Results obtained from model M5 for a t/ = 
2700 yr (see the text and table 1). A description of the 
four frames is given in the caption of Figure 1. 


that the overall morphology (particularly, the strong 
differences between the sizes and velocities of the two 
lobes) of the PPN OH231.8+4.2 can be modeled as 
the result of such an interaction. 

We have adopted a bipolar “pulse” (described by 
equations 0 El 13 and 0 instead of a continous jet, 
as an outflow of this kind appears to be implied by 
the observed Hubble-like velocity law (see §1). Also, 
we have taken the parameters (pc? and Oc)^ of 
the dense, high-velocity jet (see §2.3) that would be 
necessary for explaining the CO observations of OH 
231.8-1-4.2 of Alcolea et al. (2001). 

The parameters describing the wind from the 
companion star (i.e., QX Pup in the case of OH 
231.8+4.2) are the wind velocity {v^) and the mass 
loss rate {M^)- We have adopted a value of 40 km 
s“^ for the wind velocity, which corresponds to the 
escape velocity from QX Pup (Alcolea 2004) and 


mass loss rate values characteristic of Mira variables 
(with values ranging from 10“^ to 10 “^ Mq yr“^; see 
Table 1). We have considered a stellar wind with a 
latitude dependent density, with a density enhance¬ 
ment at the equator. The equator is perpendicular 
to the “bipolar pulse”, as suggested by the spatial 
distribution of the SiO maser associated with QX 
Pup (Sanchez Contreras et al. 2002). 

With these parameters we have computed several 
numerical simulations, from which we have obtained 
the Ha emission maps, that can be compared with 
the Ha images of OH 231.8+4.2 (Sanchez Contreras 
et al. 2000a; Bujarrabal et al. 2002). 

The numerical simulations reproduce qualita¬ 
tively well several characteristics of the Ha image 
of OH 231.8+4.2. The observed and predicted Ha 
maps have bipolar structures with strong lobe-to- 
lobe asymmetries. We should point out that the 
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Fig. 7. Results obtained for model M3 for a t/ = 1500 yr 
(see the text and table 1). The left and middle frames 
show the velocity field on the rr^^-plane (left) and ?/^-plane 
(middle). The velocity maps (arrows) are superposed on 
the density stratification, with linearly spaced isolines 
with steps of 5 X 10~^^ g cm“^. The right frame shows 
the axial velocity along the symmetry axis. 

hourglass morphology of OH 231.8+4.2 is not re¬ 
produced by these models. This discrepancy can be 
solved by including a dense disk (or torus) around 
the central source (as has been shown numerically 
by several authors; see, e.g., Erank & Mellema 1996). 
However, we have not included such a torus in the 
present simulations. 

All of our models show an upper lobe which is 
narrower and shorter than the bottom lobe, but the 
large asymmetry between the Northern and South¬ 
ern lobes of OH 231.8+4.2 is only reproduced by the 
models with a latitude dependent wind (M3 - M5). 
Eigures El island ini illustrate that the upper lobe is 
strongly reduced in size in the latitude-dependent 
wind models, and that the ratio of ^ 2 between the 
sizes of the two lobes is in qualitative agreement with 
the observations of OH 231.8+4.2 (Bujarrabal et al. 
2002). However, the width of the Southern lobe is 
better predicted by the models with spherically sym¬ 
metric winds (models Ml and M2). 

To illustrate the predicted velocity field, we 
present the velocity field of model M3 on the xz- 
plane and ^ 2 ;-plane (see EigureQ). Eor a comparison 
with the observed properties we also present the ax¬ 
ial velocity along the symmetry axis (see Eigur^ZJ- 
The predicted kinematics reproduces the strong ve¬ 
locity gradient observed along the nebular axis (Al- 
colea et al. 2001; Sanchez Contreras et al. 2000). 
However, we find that the model prediction gives ve¬ 
locities lower than the observed ones (by a factor of 
3 to 4). In the top lobe, the model produces a too 
sharp drop in the velocity. In order to obtain a bet¬ 
ter agreement with the observations, we would need 
to have a model with a higher velocity bipolar pulse. 

We have presented a number of 3D simulations 


of a bipolar pulse interacting with a continues wind 
ejected by a close companion star. In our simula¬ 
tions, we have used a dense, AGB wind with mass 
loss rates typical of AGB stars. We have shown that 
the strong asymmetries between both lobes of this 
nebula can indeed be explained in terms of a the 
proposed bipolar pulse/continous wind interaction 
model (the companion star ejecting the wind has 
to be located within the Northern lobe). We find 
that the observed lobe-to-lobe asymmetry can be ex¬ 
plained without invoking an intrinsic asymmetry in 
the ejection of the jet and the counterjet. 
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